library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggpubr)
library(rstatix)
## 
## Attaching package: 'rstatix'
## 
## The following object is masked from 'package:stats':
## 
##     filter
library(nparLD)
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:rstatix':
## 
##     select
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(PMCMRplus)
library(WRS2)
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
import seaborn as sns
import numpy as np
from scipy.stats import pearsonr, spearmanr
import os
#import decoupler as dc
import gseapy as gp
from gseapy.parser import gsea_gmt_parser
from scipy.stats import wilcoxon, ttest_ind, ttest_rel
from statsmodels.stats.multitest import multipletests
from glob import glob
from statannotations.Annotator import Annotator
plt.rcParams.update({'font.size': 14})
annot = pd.read_csv('../data/batches12_chemo_regimen_annotation.tsv', sep='\t', index_col=0)
annot['Primary Study ID'] = annot.loc[:, 'Primary Study ID'].astype(str)
pt_reg = annot.set_index('Primary Study ID').loc[:, 'chemo_regimen'].to_dict()
platinum_type_dict = annot.set_index('Primary Study ID').loc[:, 'platinum_type'].to_dict()
#tur_cys_palette = {'TUR': 'darkorange', 'Cystectomy': 'deepskyblue'}
#tur_cys_palette = {'Cystectomy': 'cornflowerblue', 'TUR': 'lightcoral'}
def dyn_comp(u, ax, lin_dens3, log_dens3, title, xlabel, test='t-test_paired', test_name = 't-test paired', size=3, format_simple = True):
  sns.boxplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    color='black', size=size)
      
  sns.lineplot(y=u, x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.35, color = 'grey')
  ax.set_yscale("log")
  if format_simple:
    ax.yaxis.set_major_formatter(matplotlib.ticker.ScalarFormatter())
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [('TUR', 'Cystectomy')], data = log_dens3, x='sample_type', y=u, order=['TUR', 'Cystectomy'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, pv)])
  annotator.annotate()
  
def dyn_comp_lin(u, ax, lin_dens3, title, xlabel, test='t-test_paired', test_name = 't-test paired', size=3):
  sns.boxplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x='sample_type', data = lin_dens3, ax=ax, order = ['TUR', 'Cystectomy'],
                    color='black', size=size)
      
  sns.lineplot(y=u, x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.35, color = 'grey')
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [('TUR', 'Cystectomy')], data = lin_dens3, x='sample_type', y=u, order=['TUR', 'Cystectomy'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, pv)])
  annotator.annotate()
  
def dyn_comp_delta(u, ax, lin_dens3, title, xlabel, to_compare, order, test = 't-test_paired', size=3, ofs = 1):
  sns.boxplot(y=u, x=to_compare, data = lin_dens3, ax=ax, order = order,
                    color = 'mediumseagreen', showfliers=False)
  sns.swarmplot(y=u, x=to_compare, data = lin_dens3, ax=ax, order = order,
                    color='black', size=size)
      
  sns.lineplot(y=u, x=to_compare, data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.7)
      
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [order], data = lin_dens3, x=to_compare, y=u, order=order)
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH', line_offset = ofs, line_offset_to_group = ofs)
  _, test_results = annotator.apply_test()._get_output()
  pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test, pv)])
  annotator.annotate(line_offset = ofs, line_offset_to_group = ofs)

Density simple comparison

norm_log_dens = pd.read_csv('../data/combined_norm_log10_density.tsv', sep='\t', index_col=0)
norm_log_dens_tumor = pd.read_csv('../data/tumor_norm_log10_density.tsv', sep='\t', index_col=0)
norm_log_dens_stroma = pd.read_csv('../data/stroma_norm_log10_density.tsv', sep='\t', index_col=0)
raw_dens = pd.read_csv('../data/all_densities_combined_kde_mm2.tsv', sep='\t', index_col=0)
dens_to_comp = ['B_cells', 'CD4_Tcells', 'Macrophages', 'Tregs', 'CD8_Tcells', 'PanCK+_cells', 'Negative_cells']
def transform_dens(norm_log_dens, addition = 0):
  norm_log_dens = norm_log_dens.loc[:, dens_to_comp] + addition
  lin_dens2 = (10**norm_log_dens).assign(sample_type = norm_log_dens.index.str.split(' ').str[1]).assign(pt_number = norm_log_dens.index.str.split(' ').str[0])
  lin_dens3 = lin_dens2[lin_dens2.pt_number.isin(lin_dens2.pt_number.value_counts()[lin_dens2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
  lin_dens3 = lin_dens3.assign(chemo_regimen = lin_dens3.pt_number.map(pt_reg)).assign(platinum_type = lin_dens3.pt_number.map(platinum_type_dict))
  log_dens2 = norm_log_dens.assign(sample_type = norm_log_dens.index.str.split(' ').str[1]).assign(pt_number = norm_log_dens.index.str.split(' ').str[0])
  log_dens3 = log_dens2[log_dens2.pt_number.isin(log_dens2.pt_number.value_counts()[log_dens2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
  log_dens3 = log_dens3.assign(chemo_regimen = log_dens3.pt_number.map(pt_reg)).assign(platinum_type = log_dens3.pt_number.map(platinum_type_dict))
  return lin_dens3, log_dens3
lin_dens3, log_dens3 = transform_dens(norm_log_dens, addition = 2)
raw_lin_dens3, raw_log_dens3 = transform_dens(raw_dens)
lin_dens3_tumor, log_dens3_tumor = transform_dens(norm_log_dens_tumor, addition = 2)
lin_dens3_stroma, log_dens3_stroma = transform_dens(norm_log_dens_stroma, addition = 2)

PDL1 simple comparison

pdl1 = pd.read_csv('../data/Prepost NAC patients overview PD-L1.csv', sep=';', index_col=0)
pdl1_cys = pdl1.loc[~pdl1.index.isna(), ['Material.1', 'CPS (%).1', 'IC (%).1', 'TC (%).1']].dropna()
pdl1_tur = pdl1.loc[~pdl1.index.isna(), ['Material', 'CPS (%)', 'IC (%)', 'TC (%)']].dropna().replace('76,5', '76.5')
pdl1_tur.index = pdl1_tur.index.astype(int).astype(str) + ' TUR'
pdl1_cys.index = pdl1_cys.index.astype(int).astype(str) + ' Cystectomy'
pdl1_cys.columns = ['Material', 'CPS (%)', 'IC (%)', 'TC (%)']
pdl1_tur = pdl1_tur.drop('Material', axis=1).astype(float)+1
pdl1_cys = pdl1_cys.drop('Material', axis=1).astype(float)+1
pdl1_df = pd.concat([pdl1_tur.assign(sample_type = 'TUR'), pdl1_cys.assign(sample_type = 'Cystectomy')])
pdl1_df_log = pd.concat([np.log10(pdl1_tur).assign(sample_type = 'TUR'), np.log10(pdl1_cys).assign(sample_type = 'Cystectomy')])
pdl1_df2 = pdl1_df.assign(pt_number = pdl1_df.index.str.split(' ').str[0])
pdl1_df3 = pdl1_df2[pdl1_df2.pt_number.isin(pdl1_df2.pt_number.value_counts()[pdl1_df2.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
pdl1_df3 = pdl1_df3.assign(platinum_type = pdl1_df3.pt_number.map(platinum_type_dict)).assign(chemo_regimen = pdl1_df3.pt_number.map(pt_reg))

pdl1_df2_log = pdl1_df_log.assign(pt_number = pdl1_df_log.index.str.split(' ').str[0])
pdl1_df3_log = pdl1_df2_log[pdl1_df2_log.pt_number.isin(pdl1_df2_log.pt_number.value_counts()[pdl1_df2_log.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False)
pdl1_df3_log = pdl1_df3_log.assign(platinum_type = pdl1_df3_log.pt_number.map(platinum_type_dict)).assign(chemo_regimen = pdl1_df3_log.pt_number.map(pt_reg))

RNA signatures simple comparison

ssgsea = pd.read_csv('../data/ssgsea_scores_Fig2.tsv', sep='\t', index_col=0)
ssgsea = ssgsea.reset_index().assign(sample_type = ssgsea.reset_index().loc[:, 'index'].apply(lambda x: x.split(' ')[1])).assign(pt_number = ssgsea.reset_index().loc[:, 'index'].apply(lambda x: x.split(' ')[0])).set_index('index')
ps_ssgsea0 = ssgsea[ssgsea.pt_number.isin(ssgsea.pt_number.value_counts()[ssgsea.pt_number.value_counts() == 2].index)].sort_values(by = ['pt_number', 'sample_type'], ascending=False)
ps_ssgsea0 = ps_ssgsea0.assign(chemo_regimen = ps_ssgsea0.pt_number.map(pt_reg)).assign(platinum_type = ps_ssgsea0.pt_number.map(platinum_type_dict))
tur_cys_palette = {'TUR': 'lightcoral', 'Cystectomy': 'lightskyblue'}
def to_cisplatin(df):
  df = df[df.platinum_type == 'Cisplatin']
  return df
lin_dens3 = to_cisplatin(lin_dens3)
log_dens3 = to_cisplatin(log_dens3)
lin_dens3_tumor = to_cisplatin(lin_dens3_tumor)
log_dens3_tumor = to_cisplatin(log_dens3_tumor)
lin_dens3_stroma = to_cisplatin(lin_dens3_stroma)
log_dens3_stroma = to_cisplatin(log_dens3_stroma)
pdl1_df3 = to_cisplatin(pdl1_df3)
pdl1_df3_log = to_cisplatin(pdl1_df3_log)
ps_ssgsea0 = to_cisplatin(ps_ssgsea0)
plt.close()
fig, axs = plt.subplots(3, 3, figsize = (15, 15))
af = axs.flat
dyn_comp('CD8_Tcells', af[2], lin_dens3,  log_dens3, '', 'CD8 T cell percentage in tumor + stroma, %')
## TUR vs. Cystectomy: t-test paired p=0.0017
dyn_comp('CD8_Tcells', af[0], lin_dens3_tumor,  log_dens3_tumor, '', 'CD8 T cell percentage in tumor, %')
## TUR vs. Cystectomy: t-test paired p=0.015
dyn_comp('CD8_Tcells', af[1], lin_dens3_stroma,  log_dens3_stroma, '', 'CD8 T cell percentage in stroma, %')
## TUR vs. Cystectomy: t-test paired p=0.0091
af[0].set_ylim((0.01), 200)
## (0.01, 200)
af[1].set_ylim((0.01), 200)
## (0.01, 200)
af[2].set_ylim((0.01), 200)
## (0.01, 200)
for i in range(3):
  af[i].set_yticks([0.01,0.1,1, 10, 100], ['0.01','0.1','1', '10', '100'])

for i,u in enumerate(['IC (%)', 'TC (%)', 'CPS (%)']):
  dyn_comp(u, af[i+3], pdl1_df3, pdl1_df3_log, '', '{} PDL1 score'.format(u), test='Wilcoxon', test_name='Wilcoxon')
  af[i+3].set_ylim((0.5, 300))
  af[i+3].set_yticks([1,10,100], ['1','10', '100'])
  

ps_ssgsea = ps_ssgsea0.copy()
dyn_comp_lin('TGFB_Mariathasan', af[6], ps_ssgsea, '', 'fibroblast TGFb ssGSEA score')
## TUR vs. Cystectomy: t-test paired p=6.5e-06
dyn_comp_lin('IFNG_signature', af[7], ps_ssgsea, '', 'INFg ssGSEA score')
## TUR vs. Cystectomy: t-test paired p=0.46
dyn_comp_lin('CD8_T_EFFECTOR', af[8], ps_ssgsea, '', 'effector CD8 T cell ssGSEA score')
## TUR vs. Cystectomy: t-test paired p=0.26
for i in range(9):
  af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
plt.subplots_adjust(hspace=0.2, wspace=0.4)
plt.show()

plt.close()
fig, axs = plt.subplots(2, 6, figsize = (22.7, 11))
af = axs.flat
cell_names = {'CD8_Tcells': 'CD8 T cell', 'Macrophages': 'Macrophage', 'CD4_Tcells': 'T helper',  'Tregs': 'Treg', 'B_cells': 'B cell', 'PanCK+_cells': 'Cancer cell', 'Negative_cells': 'Negative cell'}
for i, u in enumerate(['CD8_Tcells', 'Macrophages', 'CD4_Tcells',  'Tregs', 'B_cells', 'PanCK+_cells']):
  dyn_comp(u, af[i], lin_dens3,  log_dens3, '', '{} percentage, %'.format(cell_names[u]), test='Wilcoxon', test_name='Wilcoxon')
  af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
  
for i, u in enumerate(['CD8_Tcells', 'Macrophages', 'CD4_Tcells',  'Tregs', 'B_cells', 'PanCK+_cells']):
  dyn_comp(u, af[i+6], raw_lin_dens3,  raw_log_dens3, '', '{} density, cells/mm2'.format(cell_names[u]), test='Wilcoxon', test_name='Wilcoxon', format_simple = False)
  af[i+6].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])

plt.subplots_adjust(hspace=0.2, wspace=0.47)
plt.show()

Mixed ANOVA graphs

cys_den = log_dens3[log_dens3.sample_type == 'Cystectomy'].set_index('pt_number')
tur_den = log_dens3[log_dens3.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_den.index.intersection(tur_den.index)
delta_den = cys_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_den.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_den = delta_den.assign(pt_number = delta_den.index)
delta_den = delta_den.assign(chemo_regimen = delta_den.pt_number.map(pt_reg)).assign(platinum_type = delta_den.pt_number.map(platinum_type_dict))
cys_ps= ps_ssgsea0[ps_ssgsea0.sample_type == 'Cystectomy'].set_index('pt_number')
tur_ps= ps_ssgsea0[ps_ssgsea0.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_ps.index.intersection(tur_ps.index)
delta_ps = cys_ps.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_ps.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_ps = delta_ps.assign(pt_number = delta_ps.index)
delta_ps = delta_ps.assign(chemo_regimen = delta_ps.pt_number.map(pt_reg)).assign(platinum_type = delta_ps.pt_number.map(platinum_type_dict))
all_med = pd.read_csv('../data/all_median_distances_wo_offset.tsv', sep='\t', index_col=0)
all_med_log = np.log10(all_med.drop(['sample_type', 'pt_number'], axis=1)).assign(sample_type = all_med.index.str.split(' ').str[1]).assign(pt_number = all_med.index.str.split(' ').str[0])
all_med_lin = all_med.drop(['sample_type'], axis=1).assign(sample_type = all_med.index.str.split(' ').str[1]).assign(pt_number = all_med.index.str.split(' ').str[0])
all_med_log = all_med_log[all_med_log.pt_number.isin(all_med_log.pt_number.value_counts()[all_med_log.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False).assign(chemo_regimen = all_med_log.pt_number.map(pt_reg)).assign(platinum_type =all_med_log.pt_number.map(platinum_type_dict))
all_med_lin = all_med_lin[all_med_lin.pt_number.isin(all_med_lin.pt_number.value_counts()[all_med_lin.pt_number.value_counts() == 2].index)].sort_values(['pt_number', 'sample_type'], ascending=False).assign(chemo_regimen = all_med_lin.pt_number.map(pt_reg)).assign(platinum_type = all_med_lin.pt_number.map(platinum_type_dict))
all_med_lin = to_cisplatin(all_med_lin)
all_med_log = to_cisplatin(all_med_log)
plt.close()
fig, axs = plt.subplots(1, 2, figsize = (8, 5))
af = axs.flat
dyn_comp('CD8_Tcell_to_tumor_median', af[0], all_med_lin,  all_med_log, '', 'CD8 T cell to 1-NN cancer cell\nmedian distance, um', test='Wilcoxon', test_name='Wilcoxon')
## TUR vs. Cystectomy: Wilcoxon p=2.2e-05
dyn_comp('Macrophage_to_tumor_median', af[1], all_med_lin,  all_med_log, '', 'Macrophage to 1-NN cancer cell\nmedian distance, um', test='Wilcoxon', test_name='Wilcoxon')
## TUR vs. Cystectomy: Wilcoxon p=3.1e-05
af[0].set_yticks([10, 50, 100, 500], ['10','50','100','500'])
af[1].set_yticks([10, 50, 100], ['10','50','100'])
#af[0].set_ylim((0.0001), 2)
#af[1].set_ylim((0.0001), 2)
for i in range(2):
  af[i].set_xticks(['TUR', 'Cystectomy'], ['TUR-BT', 'Cystectomy'])
#_ = fig.set_tight_layout(True)
plt.subplots_adjust(hspace=0.2, wspace=0.5)
plt.show()

tur_dist2 = (all_med_log[all_med_log.sample_type == 'TUR']).set_index('pt_number')
cys_dist2 = (all_med_log[all_med_log.sample_type == 'Cystectomy']).set_index('pt_number')
common_index0 = tur_dist2.index.intersection(cys_dist2.index)
delta_dist0 = cys_dist2.loc[common_index0, ['CD8_Tcell_to_tumor_median', 'Macrophage_to_tumor_median']] - tur_dist2.loc[common_index0, ['CD8_Tcell_to_tumor_median', 'Macrophage_to_tumor_median']]
delta_dist0 = delta_dist0.assign(pt_number = delta_dist0.index.str.split(' ').str[0])
delta_dist0 = delta_dist0.assign(chemo_regimen = delta_dist0.pt_number.map(pt_reg)).assign(platinum_type = delta_dist0.pt_number.map(platinum_type_dict)).dropna(subset = ['chemo_regimen'])
pdl1_df3_log = pdl1_df3_log.rename(columns = {'IC (%)': 'IC', 'CPS (%)': 'CPS', 'TC (%)': 'TC'}).dropna()
pdl1_df3 = pdl1_df3.rename(columns = {'IC (%)': 'IC', 'CPS (%)': 'CPS', 'TC (%)': 'TC'}).dropna()
cys_pdl = pdl1_df3_log[pdl1_df3_log.sample_type == 'Cystectomy'].set_index('pt_number')
tur_pdl = pdl1_df3_log[pdl1_df3_log.sample_type == 'TUR'].set_index('pt_number')
common_in = cys_pdl.index.intersection(tur_pdl.index)
delta_pdl = cys_pdl.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1) - tur_pdl.loc[common_in].drop(['sample_type', 'chemo_regimen', 'platinum_type'], axis=1)
delta_pdl = delta_pdl.assign(pt_number = delta_pdl.index)
delta_pdl = delta_pdl.assign(chemo_regimen = delta_pdl.pt_number.map(pt_reg)).assign(platinum_type = delta_pdl.pt_number.map(platinum_type_dict))
def dyn_comp_anova(u, ax, lin_dens3, log_dens3, title, xlabel, loc_leg, test='t-test_paired', test_name='t-test paired'):
  sns.boxplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    color='black', dodge=True, size=3)
      
  #sns.lineplot(y=u, x='chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.5, color = 'grey')
  ax.set_yscale("log")    
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [(('MVAC', 'TUR'), ('MVAC', 'Cystectomy')), (('gemcitabine_platinum', 'TUR'), ('gemcitabine_platinum', 'Cystectomy'))], data = log_dens3, x='chemo_regimen', hue = 'sample_type', y=u, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  #pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, test_results[0].data.pvalue), '{} p={:.2g}'.format(test_name, test_results[1].data.pvalue)])
  annotator.annotate()
  handles, labels = ax.get_legend_handles_labels()
  ax.legend(handles[0:2], ['TUR-BT', 'Cystectomy'], loc=loc_leg)
  
def dyn_comp_lin_anova(u, ax, lin_dens3, title, xlabel, loc_leg, test='t-test_paired', test_name='t-test paired'):
  sns.boxplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    palette = tur_cys_palette, showfliers=False)
  sns.stripplot(y=u, x = 'chemo_regimen', hue='sample_type', data = lin_dens3, ax=ax, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'],
                    color='black', dodge=True, size=3)
      
  #sns.lineplot(y=u, hue='chemo_regimen', x='sample_type', data = lin_dens3, ax=ax, units = 'pt_number', estimator=None, linewidth=0.5, color = 'grey')
  ax.set_xlabel('')
  ax.set_ylabel(xlabel)
  ax.set_title(title)
  
  annotator = Annotator(ax, pairs = [(('MVAC', 'TUR'), ('MVAC', 'Cystectomy')), (('gemcitabine_platinum', 'TUR'), ('gemcitabine_platinum', 'Cystectomy'))], data = lin_dens3, x='chemo_regimen', hue = 'sample_type', y=u, hue_order = ['TUR', 'Cystectomy'], order = ['MVAC', 'gemcitabine_platinum'])
  annotator.configure(test=test, text_format = 'simple', comparisons_correction='BH')
  _, test_results = annotator.apply_test()._get_output()
  #pv = test_results[0].data.pvalue
  annotator.set_custom_annotations(['{} p={:.2g}'.format(test_name, test_results[0].data.pvalue), '{} paired p={:.2g}'.format(test_name, test_results[1].data.pvalue)])
  annotator.annotate()
  handles, labels = ax.get_legend_handles_labels()
  ax.legend(handles[0:2], ['TUR-BT', 'Cystectomy'], loc=loc_leg)
plt.rcParams.update({'font.size': 15})
testt = 'Wilcoxon'
testt_name = 'Wilcoxon'
plt.close()
fig, axs = plt.subplots(2, 2, figsize = (14.5, 12.5), width_ratios = [10, 5])
af = axs.flat
u = 'Macrophages'
dyn_comp_anova(u, af[0], lin_dens3, log_dens3, '', 'Macrophage cell percentage', loc_leg=(0.3, 0.03), test='Wilcoxon', test_name = 'Wilcoxon')
## gemcitabine_platinum_TUR vs. gemcitabine_platinum_Cystectomy: Wilcoxon p=0.84
## MVAC_TUR vs. MVAC_Cystectomy: Wilcoxon p=0.0011
## 
## /home/m.chelushkin/.local/lib/python3.8/site-packages/seaborn/categorical.py:166: FutureWarning: Setting a gradient palette using color= is deprecated and will be removed in version 0.13. Set `palette='dark:black'` for same effect.
##   warnings.warn(msg, FutureWarning)
dyn_comp_delta(u, af[1], delta_den, '', 'Δ log10(macrophage percentage)\nupon treatment', to_compare = 'chemo_regimen', order = ('MVAC', 'gemcitabine_platinum'), test = 'Mann-Whitney', ofs = -0.08, size=5)
## MVAC vs. gemcitabine_platinum: Mann-Whitney p=0.015
af[0].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[1].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[1].set_xlim((-0.7,1.7))
## (-0.7, 1.7)
af[1].set_ylim((-1.4,2.8))
## (-1.4, 2.8)
af[0].set_yticks([0.1,1, 10, 100], ['0.1','1', '10', '100'])
u = 'TGFB_Mariathasan'
dyn_comp_lin_anova(u, af[2], ps_ssgsea0, '', 'Fibroblast TGFb ssGSEA score', loc_leg='lower right', test='Wilcoxon', test_name = 'Wilcoxon')
## gemcitabine_platinum_TUR vs. gemcitabine_platinum_Cystectomy: Wilcoxon p=8.5e-05
## MVAC_TUR vs. MVAC_Cystectomy: Wilcoxon paired p=0.11
## 
## /home/m.chelushkin/.local/lib/python3.8/site-packages/seaborn/categorical.py:166: FutureWarning: Setting a gradient palette using color= is deprecated and will be removed in version 0.13. Set `palette='dark:black'` for same effect.
##   warnings.warn(msg, FutureWarning)
dyn_comp_delta(u, af[3], delta_ps, '', 'Δ fibroblast TGFb upon treatment', to_compare = 'chemo_regimen', order = ('MVAC', 'gemcitabine_platinum'), test = 'Mann-Whitney', ofs=0.1, size=5)
## MVAC vs. gemcitabine_platinum: Mann-Whitney p=0.038
af[2].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[3].set_xticklabels(['MVAC', 'gemcitabine +\nplatinum'])
af[3].set_xlim((-0.7,1.7))
## (-0.7, 1.7)
#_ = fig.set_tight_layout(True)
plt.subplots_adjust(hspace=0.3, wspace=0.25)
plt.show()

Mixed ANOVA models

densities

lin_dens3_cr = lin_dens3[lin_dens3.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
#lin_dens3_ptt = lin_dens3[lin_dens3.platinum_type.isin(['Cisplatin', 'Carboplatin'])]
log_dens3_cr = log_dens3[log_dens3.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
log_dens3_ptt = log_dens3[log_dens3.platinum_type.isin(['Cisplatin', 'Carboplatin'])]

CD8 T cells

py$log_dens3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(CD8_Tcells)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable   statistic     p
##   <chr>       <chr>                <chr>          <dbl> <dbl>
## 1 Cystectomy  MVAC                 CD8_Tcells     0.985 0.936
## 2 Cystectomy  gemcitabine_platinum CD8_Tcells     0.979 0.680
## 3 TUR         MVAC                 CD8_Tcells     0.969 0.499
## 4 TUR         gemcitabine_platinum CD8_Tcells     0.981 0.772
py$log_dens3_cr %>%
  group_by(sample_type) %>%
  levene_test(CD8_Tcells ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    66     0.591 0.445
## 2 TUR             1    66     1.20  0.277
box_m(py$log_dens3_cr[, "CD8_Tcells", drop = FALSE], py$log_dens3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1     0.760   0.383         1 Box's M-test for Homogeneity of Covariance Matric…
res.aov <- anova_test(
  data = py$log_dens3_cr, dv = CD8_Tcells, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F     p p<.05   ges
## 1             chemo_regimen   1  66  0.410 0.524       0.004
## 2               sample_type   1  66 11.124 0.001     * 0.048
## 3 chemo_regimen:sample_type   1  66  0.973 0.328       0.004
bwtrim(CD8_Tcells ~ chemo_regimen*sample_type, id = pt_number, data = py$log_dens3_cr)
## Call:
## bwtrim(formula = CD8_Tcells ~ chemo_regimen * sample_type, id = pt_number, 
##     data = py$log_dens3_cr)
## 
##                            value df1     df2 p.value
## chemo_regimen             0.0660   1 43.9789  0.7984
## sample_type               6.0601   1 41.1155  0.0181
## chemo_regimen:sample_type 0.4322   1 41.1155  0.5146

Macrophages

py$log_dens3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(Macrophages)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable    statistic        p
##   <chr>       <chr>                <chr>           <dbl>    <dbl>
## 1 Cystectomy  MVAC                 Macrophages     0.969 0.499   
## 2 Cystectomy  gemcitabine_platinum Macrophages     0.850 0.000156
## 3 TUR         MVAC                 Macrophages     0.881 0.00254 
## 4 TUR         gemcitabine_platinum Macrophages     0.935 0.0329
py$log_dens3_cr %>%
  group_by(sample_type) %>%
  levene_test(Macrophages ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    66     1.25  0.267
## 2 TUR             1    66     0.118 0.732
box_m(py$log_dens3_cr[, "Macrophages", drop = FALSE], py$log_dens3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1    0.0576   0.810         1 Box's M-test for Homogeneity of Covariance Matric…
bwtrim(Macrophages ~ chemo_regimen*sample_type, id = pt_number, data = py$log_dens3_cr)
## Call:
## bwtrim(formula = Macrophages ~ chemo_regimen * sample_type, id = pt_number, 
##     data = py$log_dens3_cr)
## 
##                            value df1     df2 p.value
## chemo_regimen             0.9185   1 40.8416  0.3435
## sample_type               7.7773   1 43.8798  0.0078
## chemo_regimen:sample_type 3.3520   1 43.8798  0.0739
res.aov <- anova_test(
  data = py$log_dens3_cr, dv = Macrophages, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd     F     p p<.05   ges
## 1             chemo_regimen   1  66 1.277 0.263       0.013
## 2               sample_type   1  66 9.424 0.003     * 0.041
## 3 chemo_regimen:sample_type   1  66 6.955 0.010     * 0.031
f1.ld.f1(y = py$log_dens3_cr$Macrophages, time = py$log_dens3_cr$sample_type, group = py$log_dens3_cr$chemo_regimen, subject = py$log_dens3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##            Statistic df     p-value
## Group       1.114216  1 0.291167239
## Time        9.292243  1 0.002301263
## Group:Time  7.422676  1 0.006440696

TGFb signature

ps_ssgsea_cr = ps_ssgsea0[ps_ssgsea0.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]#.drop(['15 TUR', '15 Cystectomy', '7 TUR', '7 Cystectomy', '31 TUR', '31 Cystectomy'])
py$ps_ssgsea_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(TGFB_Mariathasan)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable         statistic     p
##   <chr>       <chr>                <chr>                <dbl> <dbl>
## 1 Cystectomy  MVAC                 TGFB_Mariathasan     0.936 0.164
## 2 Cystectomy  gemcitabine_platinum TGFB_Mariathasan     0.968 0.391
## 3 TUR         MVAC                 TGFB_Mariathasan     0.926 0.102
## 4 TUR         gemcitabine_platinum TGFB_Mariathasan     0.984 0.873
py$ps_ssgsea_cr %>%
  group_by(sample_type) %>%
  levene_test(TGFB_Mariathasan ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic      p
##   <chr>       <int> <int>     <dbl>  <dbl>
## 1 Cystectomy      1    55     0.458 0.501 
## 2 TUR             1    55     6.72  0.0122
box_m(py$ps_ssgsea_cr[, "TGFB_Mariathasan", drop = FALSE], py$ps_ssgsea_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      3.38  0.0659         1 Box's M-test for Homogeneity of Covariance Matric…
res.aov <- anova_test(
  data = py$ps_ssgsea_cr, dv = TGFB_Mariathasan, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05   ges
## 1             chemo_regimen   1  55  0.945 3.35e-01       0.010
## 2               sample_type   1  55 20.482 3.26e-05     * 0.129
## 3 chemo_regimen:sample_type   1  55  3.644 6.20e-02       0.026
tgf_df = bwtrim(TGFB_Mariathasan ~ chemo_regimen*sample_type, id = pt_number, data = py$ps_ssgsea_cr)
f1.ld.f1(y = py$ps_ssgsea_cr$TGFB_Mariathasan, time = py$ps_ssgsea_cr$sample_type, group = py$ps_ssgsea_cr$chemo_regimen, subject = py$ps_ssgsea_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##             Statistic df      p-value
## Group       0.3681926  1 5.439911e-01
## Time       21.1397173  1 4.269862e-06
## Group:Time  2.2215638  1 1.360951e-01

PDL1 score

pdl_cr = pdl1_df3_log[pdl1_df3_log.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
py$pdl_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(IC)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable statistic         p
##   <chr>       <chr>                <chr>        <dbl>     <dbl>
## 1 Cystectomy  MVAC                 IC           0.797 0.0000351
## 2 Cystectomy  gemcitabine_platinum IC           0.843 0.000366 
## 3 TUR         MVAC                 IC           0.867 0.000994 
## 4 TUR         gemcitabine_platinum IC           0.890 0.00400
py$pdl_cr %>%
  group_by(sample_type) %>%
  levene_test(IC ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    61    2.47   0.122
## 2 TUR             1    61    0.0608 0.806
box_m(py$pdl_cr[, "IC", drop = FALSE], py$pdl_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      2.10   0.147         1 Box's M-test for Homogeneity of Covariance Matric…
#bwtrim(IC ~ chemo_regimen*sample_type, id = pt_number, data = py$pdl_cr)
f1.ld.f1(y = py$pdl_cr$IC, time = py$pdl_cr$sample_type, group = py$pdl_cr$chemo_regimen, subject = py$pdl_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##            Statistic df      p-value
## Group       4.987890  1 0.0255253172
## Time       13.172406  1 0.0002841017
## Group:Time  3.543588  1 0.0597760199
res.aov <- anova_test(
  data = py$pdl_cr, dv = IC, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F        p p<.05   ges
## 1             chemo_regimen   1  61  6.933 0.011000     * 0.060
## 2               sample_type   1  61 14.429 0.000338     * 0.093
## 3 chemo_regimen:sample_type   1  61  2.730 0.104000       0.019

distances

log_dist3_cr = all_med_log[all_med_log.chemo_regimen.isin(['MVAC', 'gemcitabine_platinum'])]
py$log_dist3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(CD8_Tcell_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable                  statistic         p
##   <chr>       <chr>                <chr>                         <dbl>     <dbl>
## 1 Cystectomy  MVAC                 CD8_Tcell_to_tumor_median     0.946 0.120    
## 2 Cystectomy  gemcitabine_platinum CD8_Tcell_to_tumor_median     0.957 0.159    
## 3 TUR         MVAC                 CD8_Tcell_to_tumor_median     0.779 0.0000217
## 4 TUR         gemcitabine_platinum CD8_Tcell_to_tumor_median     0.959 0.189
py$log_dist3_cr %>%
  group_by(sample_type) %>%
  levene_test(CD8_Tcell_to_tumor_median ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    66     0.495 0.484
## 2 TUR             1    66     1.76  0.190
box_m(py$log_dist3_cr[, "CD8_Tcell_to_tumor_median", drop = FALSE], py$log_dist3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      3.22  0.0728         1 Box's M-test for Homogeneity of Covariance Matric…
f1.ld.f1(y = py$log_dist3_cr$CD8_Tcell_to_tumor_median, time = py$log_dist3_cr$sample_type, group = py$log_dist3_cr$chemo_regimen, subject = py$log_dist3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##             Statistic df      p-value
## Group       0.2179797  1 6.405836e-01
## Time       20.4518025  1 6.115193e-06
## Group:Time  0.2704191  1 6.030508e-01
res.aov <- anova_test(
  data = py$log_dist3_cr, dv = CD8_Tcell_to_tumor_median, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd         F        p p<.05      ges
## 1             chemo_regimen   1  66  0.000992 0.975000       1.05e-05
## 2               sample_type   1  66 13.070000 0.000581     * 5.60e-02
## 3 chemo_regimen:sample_type   1  66  0.060000 0.807000       2.72e-04
py$log_dist3_cr %>%
  group_by(sample_type, chemo_regimen) %>%
  shapiro_test(Macrophage_to_tumor_median)
## # A tibble: 4 × 5
##   sample_type chemo_regimen        variable                   statistic        p
##   <chr>       <chr>                <chr>                          <dbl>    <dbl>
## 1 Cystectomy  MVAC                 Macrophage_to_tumor_median     0.846  4.10e-4
## 2 Cystectomy  gemcitabine_platinum Macrophage_to_tumor_median     0.927  1.81e-2
## 3 TUR         MVAC                 Macrophage_to_tumor_median     0.804  6.12e-5
## 4 TUR         gemcitabine_platinum Macrophage_to_tumor_median     0.883  1.03e-3
py$log_dist3_cr %>%
  group_by(sample_type) %>%
  levene_test(Macrophage_to_tumor_median ~ chemo_regimen)
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `data = map(.data$data, .f, ...)`.
## Caused by warning in `leveneTest.default()`:
## ! group coerced to factor.
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## # A tibble: 2 × 5
##   sample_type   df1   df2 statistic     p
##   <chr>       <int> <int>     <dbl> <dbl>
## 1 Cystectomy      1    66    0.0739 0.787
## 2 TUR             1    66    2.15   0.147
box_m(py$log_dist3_cr[, "Macrophage_to_tumor_median", drop = FALSE], py$log_dist3_cr$chemo_regimen)
## # A tibble: 1 × 4
##   statistic p.value parameter method                                            
##       <dbl>   <dbl>     <dbl> <chr>                                             
## 1      4.41  0.0357         1 Box's M-test for Homogeneity of Covariance Matric…
f1.ld.f1(y = py$log_dist3_cr$Macrophage_to_tumor_median, time = py$log_dist3_cr$sample_type, group = py$log_dist3_cr$chemo_regimen, subject = py$log_dist3_cr$pt_number, time.order = c('TUR', 'Cystectomy'), group.order = c('MVAC', 'gemcitabine_platinum'), description=F)$ANOVA.test
##  F1 LD F1 Model 
##  ----------------------- 
##  Check that the order of the time and group levels are correct.
##  Time level:   TUR Cystectomy 
##  Group level:   MVAC gemcitabine_platinum 
##  If the order is not correct, specify the correct order in time.order or group.order.

##              Statistic df      p-value
## Group       0.74679531  1 3.874928e-01
## Time       17.15151849  1 3.451334e-05
## Group:Time  0.03287889  1 8.561122e-01
res.aov <- anova_test(
  data = py$log_dist3_cr, dv = Macrophage_to_tumor_median, wid = pt_number,
  between = chemo_regimen, within = sample_type
  )
get_anova_table(res.aov)
## ANOVA Table (type III tests)
## 
##                      Effect DFn DFd      F       p p<.05      ges
## 1             chemo_regimen   1  66  1.044 0.31100       0.009000
## 2               sample_type   1  66 14.289 0.00034     * 0.080000
## 3 chemo_regimen:sample_type   1  66  0.109 0.74300       0.000662